;-------------------------------------------------------------LICENSE--------------------------------------------------------------;
;                                                                                                                                  ;
;The MAP code is written in Fortran language for magnetohydrodynamics (MHD) calculation with the adaptive mesh refinement (AMR)    ;
;and Message Passing Interface (MPI) parallelization.                                                                              ;
;                                                                                                                                  ;
;Copyright (C) 2012                                                                                                                ;
;Ronglin Jiang                                                                                                                     ;
;rljiang@ssc.net.cn                                                                                                                ;
;585 Guoshoujing Road. Pudong, Shanghai, P.R.C. 201203                                                                             ;
;                                                                                                                                  ;
;This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License         ;
;as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.             ;
;                                                                                                                                  ;
;This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of    ;
;MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.                        ;
;                                                                                                                                  ;
;You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software     ;
;Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.                                                   ;
;                                                                                                                                  ;
;-------------------------------------------------------------LICENSE--------------------------------------------------------------;

pro map_plot_2d, VAR_NAME = var_name, MESH_FLAG = mesh_flag, MAG_FIELD_FLAG = mag_field_flag, VEL_FIELD_FLAG = vel_field_flag,     $
   OUTPUTPATH = outputpath, DATAPATH = datapath, XSIZE = xsize, YSIZE = ysize, PS = ps, RGB_TABLE = rgb_table, VAR_LV = var_lv,    $
   CHAR_COLOR = char_color, LINE_COLOR = line_color, CPU_COLOR = cpu_color, BACKGROUND = background, MAG_LV = mag_lv,              $
   PLOT_RANGE = plot_range, ARROW_LEN = arrow_len, ARROW_SCL = arrow_scl, LOG = log, COLORBAR_FORMAT = colorbar_format,            $
   MAG_DIR = mag_dir, MAG_POS_START = mag_pos_start, VAR_OUT = var_out, MESH_FILL = mesh_fill, READDATA_ONLY = readdata_only,      $
   VAR_MIN = var_min, VAR_MAX = var_max, ARROW_NUM = arrow_num, MAG_PLOT_METHOD = mag_plot_method, OVERPLOT_LV = overplot_lv,      $
   TICK_LOG = tick_log, COLORBAR_FLAG = colorbar_flag, TIME_START = time_start, TIME_END = time_end, COLOR_TYPE = color_type,      $
   MAG_FIELD_SHIFT = mag_field_shift, TICKS = ticks

   if keyword_set(VAR_NAME) then begin
      var_flag = 1
   endif else begin
      var_flag = 0
      var_name = 'ro'
   endelse

   if not keyword_set(mesh_flag) then mesh_flag = 0
   if not keyword_set(mesh_fill) then mesh_fill = 0
   if not keyword_set(mag_field_flag) then mag_field_flag = 0
   if not keyword_set(vel_field_flag) then vel_field_flag = 0
   if not keyword_set(xsize) then xsize = 1200
   if not keyword_set(ysize) then ysize = 1200
   if not keyword_set(ps) then ps = 0
   if not keyword_set(mag_lv) then mag_lv = 20
   if not keyword_set(log) then log = 0
   if not keyword_set(tick_log) then tick_log = 0
   if not keyword_set(readdata_only) then readdata_only = 0
   if n_elements(arrow_num) eq 0 then arrow_num = 20
   if n_elements(var_lv) eq 0 then var_lv = 150
   if n_elements(rgb_table) eq 0 then rgb_table = -1
   if n_elements(char_color) eq 0 then char_color = 0
   if n_elements(line_color) eq 0 then line_color = 255
   if n_elements(background) eq 0 then background = 255
   if n_elements(arrow_len) eq 0 then arrow_len = 1.0
   if n_elements(colorbar_format) eq 0 then colorbar_format = '(f4.1)'
   if n_elements(time_start) eq 0 then time_start = 0
   if n_elements(time_end) eq 0 then time_end = -1
   if n_elements(arrow_scl) eq 0 then arrow_scl = 1
   if n_elements(mag_dir) eq 0 then mag_dir = 1
   if n_elements(mag_plot_method) eq 0 then mag_plot_method = 1
   if n_elements(color_type) eq 0 then color_type = 1
   if n_elements(mag_field_shift) eq 0 then mag_field_shift = 0.0
   if n_elements(ticks) eq 0 then ticks = 6

   print
   print, '|-------------------------------------------------------------------------|'
   print, ' Variable name: ', (var_flag eq 1)? var_name : 'none'
   print, ' Plot mesh? ', (mesh_flag eq 1)? 'Yes' : 'No'
   print, ' Plot magnetic field? ', (mag_field_flag eq 1)? 'Yes' : 'No'
   print, ' Plot velocity field? ', (vel_field_flag eq 1)? 'Yes' : 'No'
   print, '|-------------------------------------------------------------------------|'
   print

   usr_var, xsize, ysize
   save_info, entry_info
   usr_font

   if rgb_table eq -1 then begin
      usr_color, color = color_type
   endif else begin
      device, decomposed = 0
      loadct, rgb_table
   endelse

   cd, datapath, current = dir_old

   on_error, 2

   if file_test(datapath+'/tre.par.dac') eq 0 then begin
      print
      print, ' Error:'
      print, ' Can not find file: tre.par.dac'
      print
      retall
   endif

;----------------------------------------------------------------------------------------------------------------------------------|
   close, /all

   directory = '' & n_dims = 0 & sys = 0 & ng = 0 & ini_block_x = 0 & ini_block_y = 0 & cpu_number = 0 & output_real8_flag = 0
   nx_no_bnd = 0 & nx_block = 0 & nx_total = 0 & ny_no_bnd = 0 & ny_block = 0 & ny_total = 0 & ratio = 0 & n_levels = 0
   div_clean_flag = 0

   openr, 1, datapath + '/tre.par.dac'
   readf, 1, directory
   readf, 1, cpu_number
   readf, 1, n_dims
   readf, 1, sys
   readf, 1, ng
   readf, 1, ini_block_x
   readf, 1, ini_block_y
   readf, 1, nx_no_bnd
   readf, 1, nx_block
   readf, 1, nx_total
   readf, 1, ny_no_bnd
   readf, 1, ny_block
   readf, 1, ny_total
   readf, 1, ratio
   readf, 1, n_levels
   readf, 1, div_clean_flag
   readf, 1, output_real8_flag
   readf, 1, gm
   readf, 1, rangex
   readf, 1, rangey
   readf, 1, zero_posx
   readf, 1, zero_posy
   readf, 1, cfl
   readf, 1, min_ro
   readf, 1, max_ro
   readf, 1, min_mx
   readf, 1, max_mx
   readf, 1, min_my
   readf, 1, max_my
   readf, 1, min_mz
   readf, 1, max_mz
   readf, 1, min_bx
   readf, 1, max_bx
   readf, 1, min_by
   readf, 1, max_by
   readf, 1, min_bz
   readf, 1, max_bz
   readf, 1, min_en
   readf, 1, max_en
   readf, 1, min_az
   readf, 1, max_az
   readf, 1, min_v2
   readf, 1, max_v2
   readf, 1, min_b2
   readf, 1, max_b2
   readf, 1, min_pr
   readf, 1, max_pr
   readf, 1, min_te
   readf, 1, max_te
   readf, 1, min_beta
   readf, 1, max_beta
   close, 1

;----------------------------------------------------------------------------------------------------------------------------------|
   endian = 0
   catch, endian

   if endian ne 0 then begin
      print, '|========Old data========|'
      endian = 1
   endif

;----------------------------------------------------------------------------------------------------------------------------------|
   nt = 0l
   openr, 1, 'tre.t.dac', /f77_unformatted, swap_endian = endian
   readu, 1, nt
   t = dblarr(nt, /nozero)
   dt = dblarr(nt, /nozero)
   readu, 1, t
   readu, 1, dt
   close, 1

   catch, /cancel

;----------------------------------------------------------------------------------------------------------------------------------|
   openr, unit_block, 'tre.block.dac', /f77_unformatted, /get_lun, swap_endian = endian
   openr, unit_x, 'tre.x.dac', /f77_unformatted, /get_lun, swap_endian = endian
   openr, unit_y, 'tre.y.dac', /f77_unformatted, /get_lun, swap_endian = endian
   openr, unit_ro, 'tre.ro.dac', /f77_unformatted, /get_lun, swap_endian = endian
   openr, unit_mx, 'tre.mx.dac', /f77_unformatted, /get_lun, swap_endian = endian
   openr, unit_my, 'tre.my.dac', /f77_unformatted, /get_lun, swap_endian = endian
   openr, unit_mz, 'tre.mz.dac', /f77_unformatted, /get_lun, swap_endian = endian
   openr, unit_bx, 'tre.bx.dac', /f77_unformatted, /get_lun, swap_endian = endian
   openr, unit_by, 'tre.by.dac', /f77_unformatted, /get_lun, swap_endian = endian
   openr, unit_bz, 'tre.bz.dac', /f77_unformatted, /get_lun, swap_endian = endian
   openr, unit_en, 'tre.en.dac', /f77_unformatted, /get_lun, swap_endian = endian
   openr, unit_az, 'tre.az.dac', /f77_unformatted, /get_lun, swap_endian = endian

   n_blocks = 0l
   dx = 0.0d0
   dy = 0.0d0
   level = 0l
   nc = 0l
   block_skip = 0ll
   x_skip = 0ll
   y_skip = 0ll
   ro_skip = 0ll
   mx_skip = 0ll
   my_skip = 0ll
   mz_skip = 0ll
   bx_skip = 0ll
   by_skip = 0ll
   bz_skip = 0ll
   en_skip = 0ll
   az_skip = 0ll

   nx_base = nx_total - 2 * ng + 1
   ny_base = ny_total - 2 * ng + 1

   x = dblarr (nx_no_bnd + 1, /nozero)
   y = dblarr (ny_no_bnd + 1, /nozero)
   x_base = dblarr (nx_base, /nozero)
   y_base = dblarr (ny_base, /nozero)

   if output_real8_flag eq 8 then begin

      ro = dblarr (nx_no_bnd + 1, ny_no_bnd + 1, /nozero)
      mx = dblarr (nx_no_bnd + 1, ny_no_bnd + 1, /nozero)
      my = dblarr (nx_no_bnd + 1, ny_no_bnd + 1, /nozero)
      mz = dblarr (nx_no_bnd + 1, ny_no_bnd + 1, /nozero)
      bx = dblarr (nx_no_bnd + 1, ny_no_bnd + 1, /nozero)
      by = dblarr (nx_no_bnd + 1, ny_no_bnd + 1, /nozero)
      bz = dblarr (nx_no_bnd + 1, ny_no_bnd + 1, /nozero)
      en = dblarr (nx_no_bnd + 1, ny_no_bnd + 1, /nozero)
      az = dblarr (nx_no_bnd + 1, ny_no_bnd + 1, /nozero)

      ro_base = dblarr (nx_base, ny_base, /nozero)
      mx_base = dblarr (nx_base, ny_base, /nozero)
      my_base = dblarr (nx_base, ny_base, /nozero)
      mz_base = dblarr (nx_base, ny_base, /nozero)
      bx_base = dblarr (nx_base, ny_base, /nozero)
      by_base = dblarr (nx_base, ny_base, /nozero)
      bz_base = dblarr (nx_base, ny_base, /nozero)
      en_base = dblarr (nx_base, ny_base, /nozero)
      az_base = dblarr (nx_base, ny_base, /nozero)

      var_base = dblarr (nx_base, ny_base)

   endif else begin

      ro = fltarr (nx_no_bnd + 1, ny_no_bnd + 1, /nozero)
      mx = fltarr (nx_no_bnd + 1, ny_no_bnd + 1, /nozero)
      my = fltarr (nx_no_bnd + 1, ny_no_bnd + 1, /nozero)
      mz = fltarr (nx_no_bnd + 1, ny_no_bnd + 1, /nozero)
      bx = fltarr (nx_no_bnd + 1, ny_no_bnd + 1, /nozero)
      by = fltarr (nx_no_bnd + 1, ny_no_bnd + 1, /nozero)
      bz = fltarr (nx_no_bnd + 1, ny_no_bnd + 1, /nozero)
      en = fltarr (nx_no_bnd + 1, ny_no_bnd + 1, /nozero)
      az = fltarr (nx_no_bnd + 1, ny_no_bnd + 1, /nozero)

      ro_base = fltarr (nx_base, ny_base, /nozero)
      mx_base = fltarr (nx_base, ny_base, /nozero)
      my_base = fltarr (nx_base, ny_base, /nozero)
      mz_base = fltarr (nx_base, ny_base, /nozero)
      bx_base = fltarr (nx_base, ny_base, /nozero)
      by_base = fltarr (nx_base, ny_base, /nozero)
      bz_base = fltarr (nx_base, ny_base, /nozero)
      en_base = fltarr (nx_base, ny_base, /nozero)
      az_base = fltarr (nx_base, ny_base, /nozero)

      var_base = fltarr (nx_base, ny_base)

   endelse

   if n_elements(outputpath) then cd, outputpath

   if time_start gt nt - 1 then begin
      print, ' Warning: The given time_start is too large, which will be the last time.'
      time_start = nt - 1
   endif

   if time_end gt nt - 1 then begin
      print, ' Warning: The given time_end is too large, which will be the last time.'
      time_end = nt - 1
   endif

   if time_end eq -1 then time_end = nt - 1

;----------------------------------------------------------------------------------------------------------------------------------|
   for ti = 0, nt - 1 do begin

      if ti gt 0 then begin
         for cpui = 0, cpu_number - 1 do begin

            point_lun, unit_block, block_skip
            readu, unit_block, n_blocks

            block_skip = block_skip + (4ll + 4ll + 8ll + 8ll + 8ll) * n_blocks + 4ll + 8ll
            x_skip = x_skip + (8ll * (nx_no_bnd + 1ll) + 8ll) * n_blocks
            y_skip = y_skip + (8ll * (ny_no_bnd + 1ll) + 8ll) * n_blocks

            if output_real8_flag eq 8 then begin

               ro_skip = ro_skip + (8ll * (nx_no_bnd + 1ll) * (ny_no_bnd + 1ll) + 8ll) * n_blocks
               mx_skip = mx_skip + (8ll * (nx_no_bnd + 1ll) * (ny_no_bnd + 1ll) + 8ll) * n_blocks
               my_skip = my_skip + (8ll * (nx_no_bnd + 1ll) * (ny_no_bnd + 1ll) + 8ll) * n_blocks
               mz_skip = mz_skip + (8ll * (nx_no_bnd + 1ll) * (ny_no_bnd + 1ll) + 8ll) * n_blocks
               bx_skip = bx_skip + (8ll * (nx_no_bnd + 1ll) * (ny_no_bnd + 1ll) + 8ll) * n_blocks
               by_skip = by_skip + (8ll * (nx_no_bnd + 1ll) * (ny_no_bnd + 1ll) + 8ll) * n_blocks
               bz_skip = bz_skip + (8ll * (nx_no_bnd + 1ll) * (ny_no_bnd + 1ll) + 8ll) * n_blocks
               en_skip = en_skip + (8ll * (nx_no_bnd + 1ll) * (ny_no_bnd + 1ll) + 8ll) * n_blocks
               az_skip = az_skip + (8ll * (nx_no_bnd + 1ll) * (ny_no_bnd + 1ll) + 8ll) * n_blocks

            endif else begin

               ro_skip = ro_skip + (4ll * (nx_no_bnd + 1ll) * (ny_no_bnd + 1ll) + 8ll) * n_blocks
               mx_skip = mx_skip + (4ll * (nx_no_bnd + 1ll) * (ny_no_bnd + 1ll) + 8ll) * n_blocks
               my_skip = my_skip + (4ll * (nx_no_bnd + 1ll) * (ny_no_bnd + 1ll) + 8ll) * n_blocks
               mz_skip = mz_skip + (4ll * (nx_no_bnd + 1ll) * (ny_no_bnd + 1ll) + 8ll) * n_blocks
               bx_skip = bx_skip + (4ll * (nx_no_bnd + 1ll) * (ny_no_bnd + 1ll) + 8ll) * n_blocks
               by_skip = by_skip + (4ll * (nx_no_bnd + 1ll) * (ny_no_bnd + 1ll) + 8ll) * n_blocks
               bz_skip = bz_skip + (4ll * (nx_no_bnd + 1ll) * (ny_no_bnd + 1ll) + 8ll) * n_blocks
               en_skip = en_skip + (4ll * (nx_no_bnd + 1ll) * (ny_no_bnd + 1ll) + 8ll) * n_blocks
               az_skip = az_skip + (4ll * (nx_no_bnd + 1ll) * (ny_no_bnd + 1ll) + 8ll) * n_blocks

            endelse

         endfor

         block_skip = block_skip + 4ll + 4ll + 8ll + 8ll + 8ll

         x_skip = x_skip + 8ll * nx_base + 8ll
         y_skip = y_skip + 8ll * ny_base + 8ll

         if output_real8_flag eq 8 then begin

            ro_skip = ro_skip + 8ll * nx_base * ny_base + 8ll
            mx_skip = mx_skip + 8ll * nx_base * ny_base + 8ll
            my_skip = my_skip + 8ll * nx_base * ny_base + 8ll
            mz_skip = mz_skip + 8ll * nx_base * ny_base + 8ll
            bx_skip = bx_skip + 8ll * nx_base * ny_base + 8ll
            by_skip = by_skip + 8ll * nx_base * ny_base + 8ll
            bz_skip = bz_skip + 8ll * nx_base * ny_base + 8ll
            en_skip = en_skip + 8ll * nx_base * ny_base + 8ll
            az_skip = az_skip + 8ll * nx_base * ny_base + 8ll

         endif else begin

            ro_skip = ro_skip + 4ll * nx_base * ny_base + 8ll
            mx_skip = mx_skip + 4ll * nx_base * ny_base + 8ll
            my_skip = my_skip + 4ll * nx_base * ny_base + 8ll
            mz_skip = mz_skip + 4ll * nx_base * ny_base + 8ll
            bx_skip = bx_skip + 4ll * nx_base * ny_base + 8ll
            by_skip = by_skip + 4ll * nx_base * ny_base + 8ll
            bz_skip = bz_skip + 4ll * nx_base * ny_base + 8ll
            en_skip = en_skip + 4ll * nx_base * ny_base + 8ll
            az_skip = az_skip + 4ll * nx_base * ny_base + 8ll

         endelse

         point_lun, unit_block, block_skip
      endif

      if ti ge time_start and ti le time_end then begin

         point_lun, unit_block, block_skip
         point_lun, unit_x, x_skip
         point_lun, unit_y, y_skip
         point_lun, unit_ro, ro_skip
         point_lun, unit_mx, mx_skip
         point_lun, unit_my, my_skip
         point_lun, unit_mz, mz_skip
         point_lun, unit_bx, bx_skip
         point_lun, unit_by, by_skip
         point_lun, unit_bz, bz_skip
         point_lun, unit_en, en_skip
         point_lun, unit_az, az_skip

      endif else begin
         if ti lt time_start then continue else break
      endelse

      if ti gt time_end then break

      print, ' Time = ', t[ti]

      if readdata_only eq 0 then begin

         case var_name of
            'ro' : begin
               title = 'Density'
               range = [min_ro, max_ro]
               if n_elements(var_min) ne 0 then range[0] = var_min
               if n_elements(var_max) ne 0 then range[1] = var_max
            end
            'v2' : begin
               title = 'Total squard velocity'
               range = [min_v2, max_v2]
               if n_elements(var_min) ne 0 then range[0] = var_min
               if n_elements(var_max) ne 0 then range[1] = var_max
            end
            'b2' : begin
               title = 'Total magnetic pressure'
               range = [min_b2, max_b2]
               if n_elements(var_min) ne 0 then range[0] = var_min
               if n_elements(var_max) ne 0 then range[1] = var_max
            end
            'pr' : begin
               title = 'Gas pressure'
               range = [min_pr, max_pr]
               if n_elements(var_min) ne 0 then range[0] = var_min
               if n_elements(var_max) ne 0 then range[1] = var_max
            end
            'te' : begin
               title = 'Temperature'
               range = [min_te, max_te]
               if n_elements(var_min) ne 0 then range[0] = var_min
               if n_elements(var_max) ne 0 then range[1] = var_max
            end
            'beta' : begin
               title = 'Plasma beta'
               range = [min_beta, max_beta]
               if n_elements(var_min) ne 0 then range[0] = var_min
               if n_elements(var_max) ne 0 then range[1] = var_max
            end
            else: begin
               print, ' Wrong variable name.'
               retall
            end
         endcase

         if var_flag eq 0 then title = 'NoVar'
   
         var_levels = findgen(var_lv) / (var_lv - 1) * (range[1] - range[0] + 1d-6) + range[0]
   
         file = strcompress(title, /remove_all) + '_' +                                                                               $
               'VelF' + strcompress(vel_field_flag, /remove_all) + '_' +                                                              $
               'MagF' + strcompress(mag_field_flag, /remove_all) + '_' +                                                              $
               'Mesh' + strcompress(mesh_flag, /remove_all) + '_' +                                                                   $
               'Time' + name_seq(strcompress(t[ti], /remove_all),'3.3')
         usr_plot_sys, entry_dev, image_name = file, xsize = xsize, ysize = ysize, /start, /pixmap, ps = ps
   
         xrange = dblarr (2)
         yrange = dblarr (2)
   
         if keyword_set (plot_range) then begin
            xrange[0] = plot_range[0]
            xrange[1] = plot_range[1]
            yrange[0] = plot_range[2]
            yrange[1] = plot_range[3]
         endif else begin
            xrange[0] = 0.0d0 - zero_posx
            xrange[1] = rangex - zero_posx
            yrange[0] = 0.0d0 - zero_posy
            yrange[1] = rangey - zero_posy
         endelse
   
         contour, var_base, x_base, y_base, xstyle = 1, ystyle = 1, levels = var_levels, posi = [0.15, 0.20, 0.80, 0.85],             $
            background = background, color = char_color, /fill, title = 'Time=' + strmid(strcompress(t[ti]), 1, 4),                   $
            xtitle = 'x', ytitle = 'y', nodata = 1, xrange = xrange, yrange = yrange,                                                 $
            charsize = !usr_chars, thick = !usr_thick, charthick = !usr_chart, xthick = !usr_xthick, ythick = !usr_ythick
   
         if keyword_set(colorbar_flag) and tick_log eq 1 then colorbar, posi = [0.82, 0.20, 0.86, 0.60], /vertical,                   $
            title = title + ' (log)', /right, range = alog10(range), format = colorbar_format, background = background,               $
            color = char_color, charsize = !usr_chars_bar, charthick = !usr_chart
         if keyword_set(colorbar_flag) and tick_log eq 0 then colorbar, posi = [0.82, 0.20, 0.86, 0.60], /vertical,                   $
            title = title, /right, range = range, format = colorbar_format, background = background, color = char_color,              $ 
            ylog = log, yticks = ticks, charsize = !usr_chars_bar, charthick = !usr_chart

      endif

;----------------------------------------------------------------------------------------------------------------------------------|
      for cpui = 0, cpu_number - 1 do begin

         readu, unit_block, n_blocks

         for i_block = 0, n_blocks - 1 do begin

            readu, unit_block, level, nc, dx, dy
            readu, unit_x, x
            readu, unit_y, y
            readu, unit_ro, ro
            readu, unit_mx, mx
            readu, unit_my, my
            readu, unit_mz, mz
            readu, unit_bx, bx
            readu, unit_by, by
            readu, unit_bz, bz
            readu, unit_en, en
            readu, unit_az, az

            if nc eq 0 then begin
               case var_name of
                  'ro' : var = ro
                  'v2' : var = (mx ^ 2 + my ^ 2 + mz ^ 2) / ro ^ 2
                  'b2' : var = (bx ^ 2 + by ^ 2 + bz ^ 2)
                  'pr' : var = (gm - 1) * (en - (mx ^ 2 + my ^ 2 + mz ^ 2) / ro / 2.0 - (bx ^ 2 + by ^ 2 + bz ^ 2) / 2.0)
                  'te' : var = (gm - 1) * (en - (mx ^ 2 + my ^ 2 + mz ^ 2) / ro / 2.0 - (bx ^ 2 + by ^ 2 + bz ^ 2) / 2.0) / ro
                  'beta' : var = (gm - 1) * (en - (mx ^ 2 + my ^ 2 + mz ^ 2) / ro / 2.0 - (bx ^ 2 + by ^ 2 + bz ^ 2) / 2.0) /    $
                     (bx ^ 2 + by ^ 2 + bz ^ 2) * 2.0
                  else: begin
                     print, ' Wrong variable name.'
                     retall
                  end
               endcase

               if keyword_set(log) then begin
                  var = alog10(var)
                  range_log = alog10(range)
                  var_levels = findgen(var_lv) / (var_lv - 1) * (range_log[1] - range_log[0]) + range_log[0]
               endif

;----------------------------------------------------------------------------------------------------------------------------------|
               if readdata_only eq 0 then contour, var, x, y, levels = var_levels, /overplot, /fill, nodata = 1 - var_flag

            endif

            if readdata_only eq 0 then begin
               xb1 = min(x)
               xb2 = max(x)
               yb1 = min(y)
               yb2 = max(y)
   
               if level eq -1 then polyfill, [xb1, xb2, xb2, xb1, xb1], [yb1, yb1, yb2, yb2, yb1], color = char_color
   
               if keyword_set(cpu_color) then begin
                  if mesh_fill eq 1 then polyfill, [xb1, xb2, xb2, xb1, xb1], [yb1, yb1, yb2, yb2, yb1],                    $ 
                     color = 240 * cpui / cpu_number
                  if mesh_flag eq 1 then oplot, [xb1, xb2, xb2, xb1, xb1], [yb1, yb1, yb2, yb2, yb1],                       $
                     thick = round(!usr_thick / 2), color = 240 * cpui / cpu_number
               endif else begin
                  if mesh_fill eq 1 then polyfill, [xb1, xb2, xb2, xb1, xb1], [yb1, yb1, yb2, yb2, yb1],                    $ 
                     color = 240 * cpui / cpu_number
                  if mesh_flag eq 1 then oplot, [xb1, xb2, xb2, xb1, xb1], [yb1, yb1, yb2, yb2, yb1],                       $
                     thick = round(!usr_thick / 2), color = line_color
               endelse
            endif
         endfor
      endfor

      if readdata_only eq 0 then begin
         axis, xaxis = 0, xstyle = 1, color = char_color, xtick_get = x_tick,                                                      $
            charsize = !usr_chars, charthick = !usr_chart, xthick = !usr_xthick, ythick = !usr_ythick
         axis, yaxis = 0, ystyle = 1, color = char_color, ytick_get = y_tick,                                                      $
            charsize = !usr_chars, charthick = !usr_chart, xthick = !usr_xthick, ythick = !usr_ythick
         x_tick_len = n_elements(x_tick)
         y_tick_len = n_elements(y_tick)
         axis, xaxis = 1, xstyle = 1, color = char_color, xtickname = replicate(' ', x_tick_len),                                  $
            charsize = !usr_chars, charthick = !usr_chart, xthick = !usr_xthick, ythick = !usr_ythick
         axis, yaxis = 1, ystyle = 1, color = char_color, ytickname = replicate(' ', y_tick_len),                                  $
            charsize = !usr_chars, charthick = !usr_chart, xthick = !usr_xthick, ythick = !usr_ythick
      endif

      readu, unit_block, level, nc, dx, dy
      readu, unit_x, x_base
      readu, unit_y, y_base
      readu, unit_ro, ro_base
      readu, unit_mx, mx_base
      readu, unit_my, my_base
      readu, unit_mz, mz_base
      readu, unit_bx, bx_base
      readu, unit_by, by_base
      readu, unit_bz, bz_base
      readu, unit_en, en_base
      readu, unit_az, az_base

      pr_base = (gm - 1) * (en_base - (mx_base ^ 2 + my_base ^ 2 + mz_base ^ 2) / ro_base / 2.0 -                                  $
         (bx_base ^ 2 + by_base ^ 2 + bz_base ^ 2) / 2.0)
      te_base = pr_base / ro_base
      beta_base = pr_base / (bx_base ^ 2 + by_base ^ 2 + bz_base ^ 2) * 2.0
      jz_base = (shift(by_base, -1, 0) - shift(by_base, 1, 0)) / (x_base[2] - x_base[0]) -                                         $ 
         (shift(bx_base, 0, -1) - shift(bx_base, 0, 1)) / (y_base[2] - y_base[0])
      jz_base[0, *] = jz_base[1, *]
      jz_base[nx_base - 1, *] = jz_base[nx_base - 2, *]
      jz_base[*, 0] = jz_base[*, 1]
      jz_base[*, ny_base - 1] = jz_base[*, ny_base - 2]

      var_out = {ro:ro_base, mx:mx_base, my:my_base, mz:mz_base, bx:bx_base, by:by_base, bz:bz_base, en:en_base,                   $
         pr:pr_base, te:te_base, jz:jz_base, beta:beta_base, x:x_base, y:y_base, nx:nx_base, ny:ny_base}

      if readdata_only ne 0 then continue

;----------------------------------------------------------------------------------------------------------------------------------|
;     overplot contour
;----------------------------------------------------------------------------------------------------------------------------------|
      if n_elements(overplot_lv) ne 0 then begin
         case var_name of
            'ro' : var = ro_base
            'v2' : var = (mx_base ^ 2 + my_base ^ 2 + mz_base ^ 2) / ro_base ^ 2
            'b2' : var = (bx_base ^ 2 + by_base ^ 2 + bz_base ^ 2)
            'pr' : var = pr_base
            'te' : var = te_base
            'beta' : var = beta_base
            else: begin
               print, ' Wrong variable name.'
               retall
            end
         endcase

         annotation = replicate ('', n_elements(overplot_lv))
         for i = 0, n_elements(overplot_lv) - 1 do annotation[i] = name_seq (overplot_lv[i], '3.2')
         if keyword_set(log) then begin
            var = alog10(var)
            contour, var, x_base, y_base, levels = alog10(overplot_lv), /overplot, nodata = 1 - var_flag,                          $
               thick = round(!usr_thick / 2), c_labels = replicate(1, n_elements(overplot_lv)),                                    $
               c_charsize = !usr_chars / 2., c_charthick = !usr_chart, c_annotation = annotation
         endif else begin
            contour, var, x_base, y_base, levels = overplot_lv, /overplot, nodata = 1 - var_flag,                                  $
               thick = round(!usr_thick / 2), c_labels = replicate(1, n_elements(overplot_lv)),                                    $
               c_charsize = !usr_chars / 2., c_charthick = !usr_chart, c_annotation = annotation
         endelse
      endif

;----------------------------------------------------------------------------------------------------------------------------------|
;     magnetic field
;----------------------------------------------------------------------------------------------------------------------------------|
;      xrange = dblarr(2)
;      yrange = dblarr(2)

;      xrange[0] = min (x_base)
;      xrange[1] = max (x_base)
;      yrange[0] = min (y_base)
;      yrange[1] = max (y_base)

      if mag_field_flag eq 1 then begin
;         if mag_field_ini eq 1 and ti eq 0 then begin
;            min_az = min(az_base)
;            max_az = max(az_base)
;         endif
         if mag_plot_method eq 1 then begin
            levels = findgen(mag_lv + 1) / mag_lv * (max_az - min_az + 1d-6) + min_az + mag_field_shift

            contour, az_base, x_base, y_base, color = line_color, /overplot, levels = levels, thick = !usr_othick
         endif else begin
            vector_potential_2d, bx_base, by_base, x_base, y_base, kesai, dir = mag_dir, pos_start = mag_pos_start

            max_k = max(kesai) + 1.0d-6
            min_k = min(kesai)
            levels = findgen(mag_lv + 1) / mag_lv * (max_k - min_k) + min_k

            contour, kesai, x_base, y_base, color = line_color, /overplot, levels = levels, thick = !usr_othick
         endelse
      endif

;----------------------------------------------------------------------------------------------------------------------------------|
;     velocity field
;----------------------------------------------------------------------------------------------------------------------------------|
      if vel_field_flag eq 1 then begin
         vx_base = mx_base / ro_base
         vy_base = my_base / ro_base

         lenx = xrange[1] - xrange[0]
         leny = yrange[1] - yrange[0]
         c = 10 / leny * (rangey + rangex) / 2 / arrow_len
         devmax = convert_coord([0, 0], [0, leny / c * 1.0], /data, /to_dev)
         lenmax = sqrt((devmax[0, 1] - devmax[0, 0]) ^ 2 + (devmax[1, 1] - devmax[1, 0]) ^ 2) / 2.5
         if lenx ge leny then begin
            circle2 = fix(arrow_num)
            circle1 = fix(arrow_num) * lenx / leny
         endif
         if leny ge lenx then begin
            circle1 = fix(arrow_num)
            circle2 = fix(arrow_num) * leny / lenx
         endif
         for ii = 1, circle1 do begin
            i = fix(float(ii) * float(nx_base) / (1 + circle1)) - float(nx_base) / (1 + circle1) / 4
            for jj = 1, circle2 do begin
               j = fix(float(jj) * float(ny_base) / (1 + circle2)) - float(ny_base) / (1 + circle1) / 4

               a1 = x_base[i] + lenx / c * vx_base[i, j] / arrow_scl
               b1 = y_base[j] + leny / c * vy_base[i, j] / arrow_scl
               if x_base[i] lt xrange[0] or x_base[i] gt xrange[1] or y_base[j] lt yrange[0] or y_base[j] gt yrange[1] then continue

               devp = convert_coord([x_base[i], a1], [y_base[j], b1], /data, /to_dev)
               hsize = sqrt((devp[0, 1] - devp[0, 0]) ^ 2 + (devp[1, 1] - devp[1, 0]) ^ 2) / 2.0
               if hsize gt lenmax then hsize = lenmax

               arrow, x_base[i], y_base[j], a1, b1, color = line_color, hsize = hsize, /data, /solid, thick = !usr_othick
            endfor
         endfor

         a1 = xrange[1] + (xrange[1] - xrange[0]) / 12.0
         b1 = yrange[1] - (yrange[1] - yrange[0]) / 5.0 + leny / c * 1.0
         hsize = lenmax
         arrow, xrange[1] + (xrange[1] - xrange[0]) / 12.0, yrange[1] - (yrange[1] - yrange[0]) / 5.0,                             $
            a1, b1, color = 255 - line_color, hsize = hsize, /data, /SOLID, thick = !usr_othick
         xyouts, xrange[1] + (xrange[1] - xrange[0]) / 15.0, yrange[1] - (yrange[1] - yrange[0]) / 4.0,                            $
            '!4v!3 = ' + strcompress(arrow_scl, /remove_all), color = 255 - line_color,                                            $
            charsize = !usr_chars_bar , charthick = !usr_chart
      endif

;----------------------------------------------------------------------------------------------------------------------------------|
      usr_plot_sys, entry_dev, image_name = file, /done, ps = ps
   endfor

;----------------------------------------------------------------------------------------------------------------------------------|
   close, unit_ro
   close, unit_mx
   close, unit_my
   close, unit_mz
   close, unit_bx
   close, unit_by
   close, unit_bz
   close, unit_en
   close, unit_az

   free_lun, unit_ro
   free_lun, unit_mx
   free_lun, unit_my
   free_lun, unit_mz
   free_lun, unit_bx
   free_lun, unit_by
   free_lun, unit_bz
   free_lun, unit_en
   free_lun, unit_az
   
;----------------------------------------------------------------------------------------------------------------------------------|
   cd, dir_old
   reset_info, entry_info

end

